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Abstract 

A  lumped,  non-linear  control-oriented  dynamic  model  for  the  solid  oxide  fuel  cell  has  been  developed.  The  exponential  decay  function  and  the 
exponential  associate  function  are  introduced  to  fit  the  distribution  characteristics  of  fuel  cell  state  variables  in  the  flow  direction  of  the  gases  in 
order  to  account  for  the  effect  of  spatial  variation  of  fuel  cell  parameters  in  the  dynamic  model.  It  is  integrated  into  the  dynamic  model  by  three 
characteristic  parameters  of  the  fitting  function,  which  are  determined  via  numerical  simulations.  A  planar  solid  oxide  fuel  cell  with  co-flow  has 
been  used  to  evaluate  the  accuracy  and  applicability  of  the  current  dynamic  model.  The  dynamic  model  is  programmed  and  implemented  using  the 
SIMULINK  software.  The  simulation  results  indicate  the  model  has  good  service  quality  to  predict  the  state  variables  and  the  performance  of  the 
solid  oxide  fuel  cell. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  solid  oxide  fuel  cell  (SOFC)  is  a  promising  energy  con¬ 
version  device  because  of  its  high  efficiency,  environmental 
compatibility  and  modularity.  SOFCs  will  play  a  significant  role 
in  the  power  generation  market  in  the  future.  The  control  tech¬ 
nology  is  critical  in  its  development  and  is  an  important  factor 
for  commercializing  the  SOFC.  A  reliable  and  accurate  control- 
oriented  model  is  of  great  importance  to  understand  the  dynamic 
characteristics  of  the  SOFC. 

Control  techniques  based  on  process  models  require  simple 
and  accurate  differential  and  algebraic  equations  to  represent 
the  device  dynamic  response  characteristics.  Up  to  now,  there 
have  been  many  publications  [1-5]  on  the  dynamic  models  of 
the  SOFC.  The  characteristics  of  the  fuel  cell  have  been  sim¬ 
ulated  by  the  dynamic  models  [1,2]  at  the  cell  level  according 
to  spatial  variations.  Discretization  of  the  computational  slice 
element  along  the  streamwise  direction  was  used  to  formulate 
each  balance  of  material  and  energy.  It  is  very  difficult  to  apply 
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these  dynamic  models  to  a  controlling  study  of  the  fuel  cell 
because  the  spatial  variations  of  fuel  cell  parameters  are  included 
in  these  dynamic  models.  Lukas  et  al.  [3,4]  and  Padulles  et  al. 
[5]  introduced  non-linear  zero-dimensional,  single  node  lumped 
dynamic  models  which  are  suitable  for  fuel  cell  control  sys¬ 
tems.  However  the  spatial  effect  on  the  fuel  cell  performance 
was  not  taken  into  account  in  their  models.  The  distributions 
of  the  gaseous  molar  fractions  and  temperature  were  assumed 
to  be  uniform  within  the  whole  stack  in  the  zero- dimensional 
dynamic  model.  In  fact,  the  distribution  of  these  variables  is 
non-uniform  for  most  SOFC  configurations.  For  example,  the 
H2  molar  fraction  increases  firstly  and  then  decreases  along  the 
fuel  flow  direction  in  the  internal  reforming  SOFC.  The  perfor¬ 
mance  of  the  SOFC  depends  strongly  on  the  distributions  of  the 
gaseous  molar  fractions  and  the  temperature.  Thus,  the  accuracy 
of  the  zero-dimensional  dynamic  model  reported  in  the  literature 
[3-5]  is  not  high  for  predicting  the  fuel  cell  performances.  Comte 
et  al.  [6]  claimed  that  the  zero-dimensional  dynamic  model  could 
lead  to  as  high  as  a  20%  error  compared  with  the  more  complex 
model  in  predicting  the  thermal  and  power  outputs  for  the  planar 
SOFC  with  cross-flow.  Thus  the  spatial  effect  on  the  fuel  cell 
performance  should  be  taken  in  account  to  improve  the  accuracy 
of  the  dynamic  model. 
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Nomenclature 

the  parameter  of  exponential  fitting  function 

Ae 

fuel  cell  active  area  (m2) 

Ar 

reforming  reaction  surface  area  (m2) 

As 

anode  section  area  (m2)  (upright  with  the  gas  flow 
direction) 

Av 

volume  section  area  (m2)  (upright  with  the  gas 
flow  direction) 

CP 

specific  heat  capacity  (J  mol-1  K) 

E 

theoretical  Nernst  potential  (V) 

Ea 

activation  energy  (J  mol-1) 

fr 

the  reforming  reaction  balance  constant 

F 

Faraday  constant  (=96485  C  mol-1) 

h 

molar  enthalpy  (J  mol-1) 

AHi 

total  change  of  enthalpy  for  the  electrochemical 
reaction 

ah2 

total  change  of  enthalpy  for  the  fuel  reforming 
reaction 

A  H3 

total  change  of  enthalpy  for  the  gas-water  shift 
reaction 

i 

current  density  (A  m-2) 

io 

exchange  current  density  (A  m-2) 

ih 

diffusion  limiting  current  density  (A  m-2) 

k 

derivative  of  fitting  function  at  the  point  (0,  y(0)) 

kr,  ks 

reforming  chemical  reaction  constant 
(=4274  mol  s- 1  m2  bar),  shift  gas-water  chemical 
reaction  constant  (=1.2  x  104  molm3  s-1) 

K 

equilibrium  constant 

L 

fuel  cell  length  (m) 

M 

molecular  mass 

h 

molar  flux  (mol  s-1) 

ne 

electrons  transferred  per  reaction 

P 

operating  pressure  (Pa) 

Po 

standard  pressure  (Pa) 

Q 

molar  chemical  reaction  heat  (J  mol-1) 

Gchem 

volumetric  chemical  reaction  produced  heat 
(J  s-1  m-3) 

Ggen 

generated  heat  (J  s-1) 

r 

volumetric  reaction  rate  (mol  s-1  m-3) 

R 

universal  gas  constant  (=8.314  J  mol-1  K-1) 

Re 

electrolyte  ohmic  resistance  (£2) 

s 

source  component 

T 

temperature  (K) 

t 

time  (s) 

1 1,  t2 

fitting  function  exponential  parameter 

u 

voltage  (V) 

V 

velocity  (ms-1) 

X 

gaseous  molar  fraction 

y 

fitting  function 

z 

dimensionless  spatial  coordinate  of  gas  flow 
direction  (z  =  (z*/L)  e  [0,  1]) 

* 

z 

spatial  coordinate  of  gas  flow  direction 

Greek  letters 

0 

number  of  gas 

A 

heat  conductivity  (W  m-1  K-1) 

pi 

fluid  viscosity  (kg  ms-1) 

0 

electric  potential  (V) 

P 

density  (kgm-3) 

O 

conductivity  (S  m-1) 

% 

number  of  chemical  reaction 

K 

correction  coefficient  for  gas-water  shift  reaction 
rate  (=1.25) 

Subscripts 

a 

anode 

c 

cathode 

ch4 

methane 

CO 

carbon  monoxide 

co2 

carbon  dioxide 

e 

electrochemical  reaction 

h2 

hydrogen 

h2o 

water  (gas) 

i 

chemical 

j 

the  jth  chemical  reaction 

n2 

nitrogen 

02 

oxygen 

ohm 

ohmic 

r 

reforming  chemical  reaction 

rad 

radiant 

s 

gas-water  shift  chemical  reaction 

Superscripts 

— 

mean  value 

in 

fuel  cell  inlet 

out 

fuel  cell  outlet 

s 

solid 

Based  on  the  single  node,  a  control-oriented  dynamic  model 
of  SOFC  stack  is  developed  in  this  study.  The  distribution  char¬ 
acteristics  of  gaseous  molar  fractions  and  temperature  are  deter¬ 
mined  using  numerical  simulation  in  order  to  lump  the  spatial 
effect  into  the  model.  The  exponential  decay  function  and  the 
exponential  association  function  are  used  to  fit  the  numerical 
simulation  data  for  these  variables  in  the  flow  direction  of  the 
gases.  A  dynamic  model  with  parameters  for  the  fitting  function 
determined  from  numerical  data  is  developed  in  current  study. 
The  spatial  effect  on  the  SOFC  performance  is  considered  in 
the  present  model  via  the  parameters  in  the  fitting  functions. 
The  present  dynamic  model  has  been  programmed  in  MAT- 
LAB/S  IMULINK.  The  reliability  and  accuracy  of  the  current 
dynamic  model  is  demonstrated  through  a  planar  SOFC  dynamic 
simulation. 

2.  Dynamic  model  development 

2.7.  Chemical  processes 

The  electrochemical  process  in  a  SOFC  is  illustrated  in  Fig.  1 . 
Fuel  and  air  flow  into  the  cell  separately.  They  diffuse  through 
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Fig.  1.  Electrochemical  processes  within  SOFC. 

the  porous  electrode  structure  to  the  interlayer  and  are  then 
adsorbed.  At  the  cathode  interlayer,  oxygen  is  reduced  by  incom¬ 
ing  electrons  to  produce  oxygen  anions.  The  oxygen  anions  are 
conducted  through  the  electrolyte  to  the  anode  interlayer  where 
they  electrochemically  combine  with  the  adsorbed  hydrogen  to 
form  water  and  release  electrons  to  the  external  circuit.  The  elec¬ 
trochemical  reactions  for  the  anode  and  cathode  are  described 
as 

anode  side  :  H2  +  O2-  —>  H2O  +  2e_ 
cathode  side  :  1  /2O2  +  2e  — >►  O2 

overall  reaction  :  1  /2O2  +  H2  — >  H2O. 

For  the  internal  reforming  of  SOFC,  the  reforming  reaction 
and  gas- shift  reaction  are  included  in  the  model: 

steam-reforming  :  CH4  +  H2O  — >►  CO  +  3H2 
shifting  :  CO  +  H20=C02  +  H2 


Exponential  Decay  function: 

(a)  y  =  y0+Ale-xft'+A2e-*/t> 

Fig.  2.  The  fitting  functions  curve  character,  (a)  Exponential  decay  function:  y  = 

e-A'/fi)  4-  A2(l  —  e~x/t2). 


Also,  the  electrochemical  oxidation  of  carbon  monoxide  at 
the  anode  is  formulated  as 

CO  +  O2-^  C02  +  2e_. 

This  reaction  speed  is  2-5  times  slower  than  that  of  hydro¬ 
gen.  The  rapid  gas-shift  reaction  becomes  the  dominant  reaction 
[7,8].  Thus,  the  oxidation  of  carbon  monoxide  is  neglected  in 
the  present  model. 

2.2.  Model  development 

The  gaseous  molar  fractions  and  the  temperature  at  the 
outlet  position  of  gas  channels  are  chosen  as  state  variables 
in  the  present  study.  The  steam  reforming  reaction  is  very 
fast  at  a  short  distance  from  the  fuel  inlet  in  the  inter¬ 
nal  reforming  design  SOFC.  Consequently,  the  reactants  of 
the  reforming  reaction  CH4  and  H2O  concentrations  rapidly 
decrease  and  the  product  H2  and  CO  concentrations  increase 
within  this  area.  Then  H2  decreases  gradually  as  a  result 
of  consumption  by  the  electrochemical  oxidation  while  flow¬ 
ing  towards  the  channel  outlet.  The  exponential  decay  func¬ 
tion  and  the  exponential  associate  function  are  applied  to 
describe  the  above  characteristics  in  this  work.  The  curves 
of  the  functions  are  plotted  in  Fig.  2.  The  exponential  decay 
function  and  exponential  associate  function  are  formulated 
by  function  (1)  and  (2),  respectively.  The  distribution  of 
the  gaseous  molar  fractions  and  temperature  can  be  fitted 
through  adjusting  the  parameters  of  the  formulation  (1)  and 
(2). 

y  =  yo  +  (Aie-z/(l  +  A2c~zlt2)  (1) 

y  =  yo  +  [Ai(l  -  e-zAi )  _|_  a2(  1  -  e"^*4)]  (2) 


Exponential  Associate  function: 

(b)  y  =  To  +  A  (1 -  e  x'h  )  +  A2(\  -  e  x,tl ) 


+  Ai  e  x/t]  +  A2  e  x!t2 .  (b)  Exponential  associate  function  y  =  yo  +  Ai(l  — 
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For  the  exponential  decay  function,  there  exists  the  following 
condition  at  the  SOFC  gases  inlet: 


y(  0)  =  yo  +  Ai  +  A2 

Similarly,  at  the  outlet: 

y(l)  =  yo  +  Aie-1/;i  +  A2e~1^ 


(3) 


(4) 


At  the  point  (0,  y(0)),  the  derivative  is  defined  as  k  and  is 
given  by 


k  =  -A_A 

t\  t2 


(5) 


Combine  the  Eqs.  (3)-(5);  the  parameters  yo,  A\,  A 2  can  be 
expressed  as  follows: 


X)  —  y(0)  +  kt\  + 
t\  —  t2 


kt\(t\  -  r2)(l  -  e  1/ri) 


A 


+ 


A 


<y(0)  -  y(  l )) 


(6) 


(7) 


(8) 


t\  kthl  -e-1/'1) 

A!  =  -kt  1  -  a(y(o)  -  xd)  -  - - 

t2  kt\t2(l  —  e_1Ai) 

A2  =  -j(y(  0)  -  XD)  +  A - - 

where 

A  =  t2(  1  -  e_1/f2)  -  ti(  1  -  e"1/r>). 


Similarly,  for  the  exponential  associate  function,  the  param¬ 
eters  yo,  A i,A2  are  given  as 


yo  =  y(0) 

ti  kth  1— e_1Ai) 

Al  =  ktx  +  ^(X0)  -  XD)  +  - - 

a2  =  -^(y(0)  -  XD)  -  - - - 

where 

A 1  A2 
*  =  —  +  — 

?1  t2 

The  averaged  value  of  y  is  calculated  by 


(9) 


(10) 


(11) 


(12) 


y  =  f  y(z)dz  =  yo  -  Aiti(e  l/h  -  l)  -  A2t2(e  1/f2  -  1) 
Jo 


for  the  exponential  decay  function, 


Substituting  Eqs.  (6)— (8)  and  Eqs.  (9)— (1 1)  into  Eqs.  (13)  and 
(14),  respectively.  The  variable  y  for  both  cases  has  the  same 
form  as 


y  —  ^iy(O)  +  +  ^3 


(15) 


where 


^1  =  1  + 


t\(  1  +  t\  e  l/tl  -  t\)  -  t2(l  +  t2  e  1/?2  -  *2) 


A 


(p2  = 


t2(l  + 12  e 


-1/^2  /V) 


*2)  —  h(  1  +  X  e 


-i/fi  _ 


X) 


A 


cp2  =  &ti(l  +  ^  e 


-l/fi  __ 


h)  + 


ktl(l-e-Vtl)(l+ti  e_1Ai  —t\) 


A 


kt\t2(l  -  e  l/h){\  +  t2  e  1//2  -  t2) 


A 


2.2.7.  Mass  balance  equations 

Each  gas  material  balance  equation  is  written  as 


PV  d Xi 


RT  d t 


1  =  nmXi( 0)  -  n0UtXi(  1)  + 


(16) 


7=1 


where  i  and  j  denote  the  ith  gas  and  the  jth  independent  reactions. 
Changing  y  into  Xi  for  Eq.  (15)  and  substituting  it  into  Eq.  (16). 
then 


cp2,iPV  dxj(l) 
RT  d  t 


=  hmXi(0)  -  houtXi(  1)  +  y^yijRj.  (17) 


7=1 


Three  chemical  reactions  are  included  in  the  presented  model. 
The  calculation  of  the  reactive  rate  is  given  as  follows: 


Ri=A{ 


nPF 


R2  =  ArkrfrP  I  A'CII4  exp  ( |  d z 


RT 


(18) 


(19) 


l 


7?3  =  $AsksP  /  (XC0*H20  -  ^shift-^C02xH2)dz 

0 


(20) 


T^shift  is  obtained  using  the  method  of  Ref.  [9].  The  outlet  of 
total  molar  number  is  calculated  by  the  continuity  equations  as 
described  in  the  following  forms: 


(13)  fiout  =  nm  +  27?2  for  the  anode  side, 


y  =  [  y(z)dz  =  yo  ~  Ai(l  +  t\  e 
Jo 

—  A2(  1  +  t2  e_1///2  —  t2) 


nout  =  hm - R\  for  the  cathode  side. 

2 

The  matrix  form  for  the  mass  balance  equations  is  written  as 
follows: 


1  . 


for  the  exponential  associate  function. 


(14)  Cx  =  Ax  +  B 


(21) 
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where 


PV 

C  = 

RT 


<P2,\ 

0 

0 

0 

0 


0 

<P2,2 

0 

0 

0 


0 

0 

<P2,3  •  •  • 

0 

0 


out 


1 

0 

0 


0  0 
1  0 
0  1 


0 

0 

0 


0  0  0  0 
0  0  0  ...  1 


0 

0 

0 

0 

^2,7] 


V  k 


Mscsp  +  (  Av^c^ 


k=  l 


i=\ 


df 
d  t 


2  rj  k 


Cpi[nfxi(0)T(0)  -  h0UtXj(l)T(l)] 


k=  1  i=  1 
2  m 


k=  1  z=l 


-  d  n; 

hiU  +  Q 


gen' 


(24) 


The  response  time  of  gas  molar  fraction  is  much  faster  than  that 
of  temperature.  Therefore,  the  second  term  on  the  right-hand 
side  of  Eq.  (24)  can  be  neglected.  Expanding  the  left  terms  and 
the  Eq.  (24)  becomes 


<P2,T 


dT(\) 
d  t 


I 

xi(0) +  '52njRj 

7=1 

I 

*2(0)  +  ^2y2,jRj 
7=1 


2  r)k 


k=  1  i= 1 


hfxi (0)7(0)  -  «r^'(l)r(l) 


(25) 


where 

2  gen  =  +  AH2R2  +  AH3R3  —  Edc- 


I 

*^(0)  +  E^7^7 
7=1 


2.2. 3.  Voltage  model 

The  Nernst  potential  is  calculated  as  a  function  of  current 
density,  operating  pressure,  temperature  and  the  gases  molar 
fraction.  The  irreversibility  in  the  voltage  drop  is  contributed  by 


T  =  cp\jT(  0)  +  ^2,r^(l)  +  cp2j  (22) 

where 

ri 

y>(D  =  1 

i=l 


2.2.2.  Energy  balance  equation 

The  temperature  dynamic  model  is  acquired  from  the  energy 
conservation  equation.  The  assumption  that  the  temperature  dis¬ 
tribution  meets  the  one-dimensional  distribution  along  the  flow 
direction  is  similar  to  the  molar  fractions  assumption.  The  gases 
thermal  capacity  is  calculated  separately  for  the  fuel  and  the  air. 
The  energy  conservation  equation  is  given  as 


s 


Me 


p 


df 

d7 


2  m  -  _ 

v-^v-^d  bin 


+  EE 

k=  1  i=  1 


rn 


2  ilk 


d  t 


=  [Ai(0)*/(0)  -  A,-(1)A/(1)] 

k=  1  i= 1 


+  Q 


gen 


(23) 


where 

771  =  5  for  H2,  H20,  CO,  C02,  CH4  while  k=  1,  r)2  =  2  for  02, 
N?  while  /:  =  2. 

Assuming  the  gas  specific  heat  is  a  function  of  the  average 
temperature  f  and  the  thermodynamic  properties  of  the  mixture 
gases  satisfy  the  ideal  gas  mixture  law.  Then 


Channel 
reaction  area 


Y  Fuel 

Symmetry  surface 


Under  rib 
reaction  area 


Interconnector 


Fig.  3.  The  geometry  of  the  half  single-unit  cell  for  the  co-flow  planar  SOFC. 
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Table  1 


The  diffuse  coefficient  and  source  terms  for  the  governing  equations 


Diffuse  coefficient  (T^) 

Source  terms  (Sa) 

Electrode 

Interlayer 

Mass 

0 

anode  side  :  Mh2oKh9o,i^i  +  47h9Kh9,i^i,  cathode  side  :  Mo2Ko9,i^i 

N-S 

~(n/a)V 

pD/,m 

Mi(yU2R2  +  Yi,2R3) 

Mi(yiARi  +  yi,2R2  +  y^Rf) 

Electric  potential 

a 

anode  side:  —i,  cathode  side:  i 

Energy 

/jl/P  r 

(2ohm  +  (2  rad 

Qo\xm  +  Qc  hem  +  (2  rad 

Table  2 

Material  properties  and  operating  conditions  [11,12] 


Cathode 

Electrolyte 

Anode 

Interconnect 

<7 

4-2y°7exp(  i2r°°) 

3.34  x  104  exp  (  10200 ) 

9-5xr107  exp  (  nr50) 

9-3*rlo6exp(  1100  ) 

7 

3 

2 

3 

3.5 

P 

6600 

cp 

400 

Operating  conditions 
Operating  voltage  (V) 

0.65 

Operating  pressure  (bar) 

1 

Inlet  gas  temperature  (K) 
Inlet  gas  composition 

Air:  02  21%,  N2  79%;  fuel:  H2  26.26%,  H20 
49.34%,  CO  2.94%,  C02  4.36%,  CH4  17.10% 

1173 

Table  3 

comparisons  of  the  numerical  simulation  results  with  IEA  benchmark 


Voltage 

Average  current  density 

Maximum  temperature 

Minimum  temperature 

Fuel  utilization 

Benchmark 

0.633-0.649 

3000 

1294-1307 

1120-1135 

85% 

Simulation 

0.65 

2959 

1301 

1139 

83% 

Table  4 

The  fitting  function  parameters  for  each  of  state  variables 


Xh2 

Xo2 

*h2o 

*co 

*co2 

*ch4 

Vn2 

T 

h 

0.12012 

0.98607 

0.10457 

0.08696 

23.20000 

0.10269 

9.60000 

0.08372 

h 

0.80000 

0.29651 

1.82261 

1.65348 

17.90000 

0.03254 

9.50000 

0.92091 

k 

2.90000 

-0.01510 

-2.30000 

1.22000 

0.11564 

-2.50000 

0.02682 

-815.00000 

three  kinds  of  polarization  losses  named  the  activation,  concen¬ 
tration  and  ohmic  overpotentials.  The  Nernst  equation  for  the 
SOFCis: 


RT 

E=—\n  K(T ) 
nQF 


2  nQF 


xpl2oP° 

xfxo2 


(26) 


where  Pa  =  Pc  =  P  is  assumed.  Subtracting  the  electrochemical 
losses  and  the  output  voltage  is  given  by 


U  =  E  —  lRQ(T) 


1R^~  sinh-1  ( — ) 
nQF  V2*0a/ 


2  RT 
nQF 


sinh 


-l 


+ 


RT 


In 


(27) 


Table  5 

Comparison  of  the  value  of  dynamic  model  prediction  with  numerical  simulation  for  the  variables 


Xh2 

Xo2 

*h2o 

*co 

•*C02 

•*CH4 

Xn2 

T 

U 

Dynamic  model  prediction 

0.08384 

0.1854 

0.7343 

0.0563 

0.1255 

6.0e— 5 

0.8146 

1304 

0.655 

Numerical  simulation  (average  values) 

0.1124 

0.183 

0.7055 

0.027 

0.1549 

2.0e— 4 

0.817 

1301 

0.65 
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where  the  second  term  on  the  right-hand  side  of  the  above  equa¬ 
tion  is  ohmic  loss,  the  third  term  is  activation  loss  at  the  anode 
side,  the  fourth  term  is  activation  loss  at  the  cathode  side,  and 
the  last  term  is  concentration  loss. 


3.  Model  evaluation 

The  dynamic  response  performance  of  a  planar  SOFC  with 
co-flow  designed  with  the  direct  internal  reforming  method  is 


Fig.  4.  Non-linear  curve  fitting  using  the  exponential  decay  function  and  exponential  association  function  for  the  gas  molar  fraction  and  the  temperature  along  the 
gases  flow  direction  in  the  gas  channel:  (a)  H2;  (b)  O2;  (c)  CO;  (d)  CH4;  (e)  H2O;  (f)  CO2;  (g)  N2;  (h)  temperature. 
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(h)  z/L 


Fig.  4.  ( Continued ). 


simulated  to  demonstrate  the  reliability  and  accuracy  of  the 
present  dynamic  model.  This  fuel  cell  stack  was  designed  as 
a  modeling  excise,  which  was  launched  by  the  International 
Energy  Agency  (IEA)  involving  seven  Europeans  countries  and 
Japan  in  1995  [12]. 

The  spatial  distributions  of  state  variables  in  the  SOFC 
are  calculated  by  solving  the  three-dimensional  Navier-Stokes 
equations  combined  with  gases  transport  equations.  The  numer¬ 
ical  results  are  used  to  find  the  parameters  t\ ,  t2 ,  and  k  of  the 
fitting  function  in  the  present  dynamic  model.  Different  values 
of  1 1,  t2 ,  and  k  are  obtained  for  each  gas  species  and  tempera¬ 
ture  spatial  distribution  in  this  case.  The  numerical  method  for 
solving  three-dimensional  Navier-Stokes  equations  and  gases 
transport  equations  is  briefly  described  in  later  section. 

3.1.  Numerical  simulations 

The  governing  equation  for  the  steady  operation  can  be  writ¬ 
ten  as  follows: 

di v(pV(p)  =  div(/0  grad 0)  +  +  Sa  (28) 

where  0  is  a  generic  variable. 

The  diffusion  coefficient  and  source  terms  Sa  are  summarized 
in  Table  1 .  The  source  term  Sa  for  energy  equation  is  the  sum 
of  the  ohmic  heat,  the  chemical  reaction  heat  and  the  radiant 
heat  exchange.  The  discrete  ordinate  radiation  model  is  applied 
to  calculate  the  radiant  heat  exchange  [10]  in  present  numeri¬ 
cal  simulation.  The  ohmic  heat  is  calculated  by  the  following 
formulae: 

(2 ohm  =  or  grad  0  grad  0  (29) 

At  the  interlayer,  the  chemical  reaction  heat  is  calculated  by 

Gchem  —  Ge  T  A  Gr  T  Fs  Gs  (20) 

The  electrochemical  thermodynamic  models  are  shown  in 
Eqs.  (26)  and  (27). 

The  given  parameters  in  present  numerical  simulation  are 
specified  in  Table  2.  The  flow  is  assumed  to  be  laminar  because 


the  Reynolds  number  is  very  low  ( Re<  100)  in  this  case.  The 
governing  equations  including  continuity,  momentum  conser¬ 
vation,  transport,  energy  conservation,  and  electric  potential  are 
solved  using  the  computational  fluid  dynamics  software  Fluent 

6.1.  Thermodynamic  equations  for  the  electrochemical  reactions 
are  incorporated  and  solved  by  user-defined  functions  (UDFs), 
which  are  provided  by  Fluent  6.1.  The  governing  equations  are 
discretized  with  second  order  upwind  format  using  control  vol¬ 
ume  method.  Half  of  the  single-unit  is  simulated  because  of 
symmetry. 

Half  of  the  single-unit  cell  planar  SOFC  is  illustrated  in 
Fig.  3.  The  geometrical  parameters  of  the  single-unit  cell  pla¬ 
nar  SOFC  are  given  in  [11]  and  are  summarized  as  follows: 
active  area:  5.42  x  100  mm,  anode  thickness:  50  pan,  cathode 
thickness:  50  p,m,  electrolyte  thickness:  150  |jim,  channel  height: 
1  mm,  rib  width:  2.42  mm,  and  bipolar  plate  thickness:  2.5  mm. 
The  permeability  of  porous  for  electrode  is  1 .0  x  10~8  m-2.  And 
its  porosity  is  0.2. 

Table  3  gives  the  comparison  of  current  numerical  simulation 
results  with  the  data  given  by  the  IEA.  As  shown  in  Table  3, 
the  numerical  simulation  results  agree  with  the  IEA  benchmark 
values.  This  indicates  that  the  accuracy  of  the  results  obtained 
by  numerical  simulation  with  current  model  is  reliable. 

3.2.  Dynamic  simulation 

The  parameters  of  fitting  functions  t\ ,  t2,  and  k  must  be  deter¬ 
mined  before  the  present  dynamic  model  can  be  implemented. 
The  kind  of  functions  used  in  the  curve  fitting  is  decided  accord¬ 
ing  to  the  distribution  characteristics  of  each  state  variable.  In 
this  paper,  only  the  state  variable  distribution  characteristics 
within  the  channel  area  are  fitted  by  the  introduced  fitting  func¬ 
tions.  The  exponential  decay  function  is  employed  for  the  molar 
fraction  of  H2,  O2,  CO,  CH4  The  exponential  association  func¬ 
tion  is  employed  for  the  molar  fractions  of  H2O,  CO2,  N2  and 
the  temperature. 

Fig.  4(a)-(d)  shows  the  fitting  curves  of  molar  fraction  of  H2, 
O2,  CO,  CH4  using  the  exponential  decay  function.  Fig.  4(e)-(h) 
shows  the  fitting  curves  of  the  molar  fractions  of  H2O,  CO2,  N2 
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and  the  temperature  using  the  exponential  association  function. 
The  solid  lines  in  Fig.  4  are  the  fitted  functions  using  the  numeri¬ 
cal  simulation  results.  As  shown  in  Fig.  4,  the  fitting  curves  well 
represent  the  numerical  results  with  the  solid  curves  very  close  to 
the  numerical  data  at  all  points,  which  represent  the  distribution 


characteristics  of  state  variables  along  the  gases  streamwise.  The 
fitting  functions  parameters  A\,  A2,  t\ ,  and  t2  for  each  of  state 
variables  are  then  determined.  The  parameter  k  is  calculated  by 
Eq.  (5)  for  molar  fraction  of  H2,  O2,  CO,  CH4  and  by  Eq.  (12) 
for  the  molar  fractions  of  H2O,  CO2,  N2  and  the  temperature, 


Fig.  5.  The  response  of  state  variables  and  the  performance  for  the  planar  co-flow  SOFC:  (a)  H2;  (b)  O2;  (c)  H2O;  (d)  CO;  (e)  CO2;  (f)  CH4;  (g)  temperature;  (h) 
voltage. 
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(g)  t  (sec) 


(h)  t (sec) 


Fig.  5.  ( Continued ). 


respectively.  The  values  of  t\ ,  U  and  k  for  each  of  state  variable 
are  presented  in  Table  4. 

The  dynamic  simulation  is  implemented  using  SIMULINK. 
The  inlet  gaseous  molar  fractions,  the  inlet  temperature,  the 
averaged  current  density,  and  the  fuel  utilization  are  specified 
as  input  parameters.  The  material  properties  are  identical  to 
the  numerical  simulation  used.  Table  5  shows  the  comparison 
of  steady  value  of  present  dynamic  simulation  with  the  data 
obtained  by  numerical  simulation.  In  Table  5,  the  numerical  data 
are  the  outlet  average  values  of  each  state  variable.  The  voltage 
is  the  fuel  cell  output  voltage.  It  can  be  seen  from  Table  5  that 
the  value  of  dynamic  model  prediction  for  xo2 ,  ^n2  >  ^ch4  ,  T 
and  U  is  consistent  with  that  of  numerical  solution.  The  dif¬ 
ferences  between  the  value  of  dynamic  model  prediction  and 
the  numerical  simulation  results  for  xn2,  yh2o,  yco  and  vco2 
are  less  than  0.03  comparing  the  numerical  simulation  results. 
This  difference  is  caused  by  the  interconnector  rib  effect.  As 
shown  in  Fig.  3,  for  the  underneath  interconnector  rib  reaction 
area,  the  gaseous  species  are  transported  only  by  diffusion  in 
the  porous  electrodes.  This  leads  to  a  big  difference  of  gaseous 
molar  fractions  between  the  underneath  interconnector  rib  reac¬ 
tion  area  and  the  channel  reaction  area.  In  the  dynamic  model, 
the  parameters  of  fitting  functions  only  reflect  the  distribution 
of  gaseous  molar  fractions  in  the  channel  reaction  area.  This 
causes  errors  in  the  prediction  of  gas- shift  reaction  rate.  Present 
dynamic  model  only  uses  a  correction  coefficient  f  in  Eq.  (20) 
to  correct  the  gas- shift  reaction  rate.  This  cannot  completely 
reflect  the  difference  of  gaseous  molar  fractions  between  the 
channel  reaction  area  and  the  underneath  interconnector  reaction 
area. 

Fig.  5  shows  the  transient  response  of  state  variables  for  the 
steady-state  operation.  Most  of  the  state  variables  are  first  order 
response  except  the  CO  molar  fraction  which  is  a  second  order 
term.  The  time  constant  of  thermal  response  is  much  larger  than 
that  of  gaseous  molar  fractions.  And  the  output  voltage  has  the 
minimum  time  constant. 

4.  Conclusion 

A  non-linear  control-oriented  dynamic  model  has  been 
developed  based  on  the  assumption  that  the  variables  follow 
a  one-dimensional  distribution.  Two  kinds  of  fitting  function, 


namely  the  exponent  decay  function  and  the  exponent  associate 
function  were  introduced  to  fit  the  distribution  characteristics 
of  the  gaseous  molar  fractions  and  temperature  along  the 
streamwise  direction.  The  spatial  effect  has  been  lumped  into 
the  dynamic  model  by  fitting  the  three  parameters  t\ ,  t2,  and  k 
of  the  used  function.  These  parameters  are  determined  through 
numerical  simulations. 

A  co-flow  configuration  planar  SOFC  has  been  introduced 
to  evaluate  the  developed  model.  The  dynamic  model  has  been 
implemented  on  this  planar  SOFC  using  SIMULINK  software. 
The  simulation  results  show  that  the  present  dynamic  model  has 
a  good  predictive  quality  for  the  state  variable  responses  and 
fuel  cell  performance. 

Acknowledgments 

The  research  work  was  sponsored  by  the  Doctor  Foundation 
of  Xi’an  Jiaotong  University  (No.  DFXJTU2005-01)  and  the 
National  High  Technology  Research  and  Development  Program 
of  China  (No.  2002AA503020). 

References 

[1]  R.S.  Gemmen,  E.  Liese,  J.G.  Rivera,  F.  Jabbari,  J.  Brouwer,  Proceedings 
of  the  ASME  Turbo  Expo,  Munich,  Germany,  2000  (GTO-554). 

[2]  R.A.  Roberts,  F.  Jabbari,  J.  Brouwer,  R.S.  Gemmen,  E.A.  Liese,  Proceed¬ 
ings  of  the  ASME  Turbo  Expo,  Georgia,  USA,  2003  (GT-38774). 

[3]  M.D.  Lukas,  K.Y.  Lee,  H.  Ghezel-Ayagh,  IEEE  Trans.  Energy  Convers.  14 
(1999)  1651. 

[4]  M.D.  Lukas,  K.Y.  Lee,  H.  Ghezel-Ayagh,  IEEE  Trans.  Energy  Convers.  16 
(2001)  289. 

[5]  J.  Padulles,  G.W.  Ault,  J.R.  McDonald,  J.  Power  Sources  86  (2000)  495. 

[6]  A.  Comte,  Ph.  Matheron,  Ph.  Stevens,  Proceedings  of  the  Third  European 
SOFC  Form,  Nantes,  France,  1998. 

[7]  R.J.  Braun,  Optimal  Design  and  Operation  of  Solid  Oxide  Fuel  Cell  Sys¬ 
tems  for  Small-scale  Stationary  Applications,  PhD  Thesis,  University  of 
Wisconsin-Madison,  USA,  2002. 

[8]  Y.  Matsuzake,  I.  Yasuda,  J.  Electrochem.  Soc.  147  (2000)  1630. 

[9]  H.  Yakabe,  T.  Ogiwara,  M.  Hishinuma,  I.  Yasuda,  J.  Power  Sources  102 
(2001)  144. 

[10]  Fluent  User’s  Guide,  Version  6.1,  Fluent  Inc.,  Lebanon,  NH,  2003. 

[11]  J.R.  Ferguson,  J.M.  Fiard,  R.  Herbin,  J.  Power  Sources  58  (1996)  106. 

[12]  E.  Achenbach,  SOFC  stack  modeling,  Final  Report  of  Activity  A2,  Annex 
II:  Modeling  and  Evaluation  of  Advanced  Solid  Oxide  Fuel  Cells,  Inter¬ 
national  Energy  Agency  Programme  on  R,  D&D  on  Advanced  Fuel  Cells, 
Juelich,  Germany,  1996. 


